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2-D/Axisymmetric Formulation of 
Multi-dimensional Upwind Scheme 


William A. Wood* and William L. Kleb* 

NASA Langley Research Center, Hampton, VA 23681 

A multi-dimensional upwind discretization of the two-dimensional/axisymmetric 
Navier-Stokes equations is detailed for unstructured meshes. The algorithm is an ex- 
tension of the fluctuation splitting scheme of Sidilkover. Boundary conditions are im- 
plemented weakly so that all nodes are updated using the base scheme, and eigen-value 
limiting is incorporated to suppress expansion shocks. Test cases for Mach numbers rang- 
ing from 0.1—17 are considered, with results compared against an unstructured upwind 
finite volume scheme. The fluctuation splitting inviscid distribution requires fewer op- 
erations than the finite volume routine, and is seen to produce less artificial dissipation, 
leading to generally improved solution accuracy. 


Nomenclature 

A Flux Jacobian 

A Auxiliary variables flux Jacobian 

B Axisymmetric source term 

Cf Skin friction coefficient 

C p Pressure coefficient 

D Linearity preserving matrix 

E Total energy 

F Flux function 

H Total Enthalpy 

M Upwinding matrix 

P Pressure 

Q Limiter ratio 

R Gas constant 

R e Reynolds’ number 

S Area 

T Temperature 

U Conserved variables 

V Primitive variables 

V Projected velocity 

W Auxiliary variables 

X Eigen-vectors 

X Auxiliary variables eigen-vectors 

Z Parameter vector 

a Sound speed 

c p Specific heat 

e Internal energy 

t Length 

n,t Normal/tangential vectors 

q Heat-transfer rate 

r Position vector 

t Time 
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u,v 

Velocity components 

x,y 

Cartesian coordinates 

r 

Generalized integration surface 

A 

Eigen-values 

n 

Generalized integration volume 

$ 

Finite volume artificial dissipation 

<*,0 

Curvilinear advection speeds 

S 

Incremental amount 

7 

Ratio of specific heats 

K 

Thermal conductivity 

p 

Coefficient of viscosity 

<t> 

Fluctuation 

V’ 

Limiter function 

p 

Density 

T 

Shear stress 

V 

Finite element shape function 

w 

Axisymmetric switch 

Cv 

Curvilinear coordinates 

Superscripts: 

i 

Inviscid 

V 

Viscous 

T 

Transpose 

x,y,£,r) 

Spatial component of a vector 

* 

Second-order fluctuation 

/ 

Fluctuation splitting artificial dissipation 

Subscripts: 

0 

Current node 

oo 

Freestream 

0 

Stagnation value 

w 

Wall 

R,L 

Right/left 

T 

Triangle 

i.i.k 

Indicies 
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Notation: 

AOA Angle of attack 

COE Contributions from other elements 
LHS Left-hand side 

RHS Right-hand side 

V Gradient 

A Backward difference 

e Permutation operator 


Bold indicates vectors of the system of equations. 
The vector symbol, indicates spatial vectors. Tilde 
quantities are Roe-averaged, while the overbar is for 
linearly averaged quantities. Hats denote unit vectors. 
The breve symbol, ", indicates quantities in auxiliary 
variables. Subscripts of other variables indicate differ- 
entiation. 


Introduction 

U PWIND fluctuation splitting and finite volume 
discretization schemes are detailed for the two- 
dimensional and axisymmetric equations of motion for 
a perfect gas on triangulated domains. Both the finite 
volume and fluctuation splitting upwind schemes are 
applied to the inviscid flux, while the viscous flux is 
discretized with a scheme analogous to finite element. 

Verification and validation of the schemes is per- 
formed using the test cases and methodology of Shing- 
hal 1 and Roache, 2 with examples ranging from the 
incompressible flat plate to a Mach- 17 cylinder. 

New contributions include the axisymmetric for- 
mulation of the upwind fluctuation splitting distribu- 
tion, the proper form for eigen-value limiting for this 
scheme, the head-to-head comparison of finite volume 
and fluctuation splitting, and the application of fluctu- 
ation splitting to a hypersonic heat-transfer validation 
test. 


Formulations 

The Navier 3 -Stokes 4 system of equations can be 
written in two-dimensional or axisymmetric non- 
dimensional form as, 


The inviscid and viscous fluxes are, 


' pV N 

puV + (1,0)P 
pvV + (0, 1)P 
V pVH ) 



°° \kVT + VtJ 

with the shear-stress tensor defined, 


r =p 


v T v + (v T v) T - VI 


( 4 ) 


( 5 ) 


(6) 


The inviscid and viscous axisymmetric source terms 
each have only one non-zero term, 

Bi=P - = <7) 

The governing equations are discretized using two 
different, second-order node-based schemes for un- 
structured (triangulated) meshes. The popular 
Barth 5,6 finite volume scheme is chosen as the baseline 
for comparison. The other scheme is the multi-dimen- 
sional-upwind fluctuation splitting discretization due 
to Sidilkover, 7 9 extended here to include eigen-value 
limiting, axisymmetric terms, and both thin-layer and 
full Navier-Stokes viscous terms. 


State Vector 

In the finite volume context, integration of the de- 
pendent variables over the control volume about node i 
is performed as, 



t dEl — Z*7 a Sj\J i t 


(8) 


For two-dimensional w a = 1, while for axisymmetric 
w a can be either taken as w a = ip, for mass-lumping 
to the node, or as the y - value of the centroid of f!., . 

In the fluctuation splitting context, the parameter 
vector is taken to vary linearly over each element. For 
a perfect gas, changes to the conserved variables can 
be related to changes in the parameter vector as, 


zu a V t +W-(zu a F i ) = V-(zu a F v )+wB i -wB v ( 1 ) 

where w is a logical switch between two-dimensional 
(w = 0) and axisymmetric (zu = 1) equations and, 

zu a = 1 — w + zuy (2) 


U z 


d U = TJz dZ 

2Z, 0 0 0 

Z 2 Z\ 0 0 

Z% 0 Z\ 0 

~z 4 ^ Z 2 Z 3 i Zy 


(9) 

(10) 


is 1 for two dimensions and y for axisymmetric. 

The conserved state vector is, 

U = (p, pu , pv, pE) T (3) 


Integration of u7 a U * over an element leads to a mass 
matrix, 


/ ZUaVtdEl = / WaBz'ZtdVl 

Jq Jq 


(ID 
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If mass-lumping to the nodes is employed, introduc- 
ing temporal, but not spatial, errors, Eqn. 11 can be 
partitioned among the three nodes defining Q as, 

[ WaXJtdn = ^ (12) 

Ja 3 i=l 

so that the sum of all contributions to node i equals 
Axisymmetric Sources 

In the finite volume framework, the inviscid axisym- 
metric source term, B\ is simply evaluated at the node 
as, 



B' dil =Sj B) 


(13) 


While some authors insist on upwinding source 
terms for fluctuation splitting , 8,10 the present analysis 
considers an upwind distribution to be inappropriate 
for the axisymmetric source terms, which arise from 
purely geometric manipulations. The axisymmetric 
source term can be distributed to the node in the fluc- 
tuation splitting framework following a mass-lumped 
analogy as, 


w aj Sj\J jt e- zuS 3 Bj (14) 


which is equivalent to the finite volume treatment of 
Eqn. 13. A modification of this distribution is to send 
contributions weighted by the averaged values, 

WajSj V it ^w^Wj + COE (15) 

A more rigorous treatment integrates the source 
term analytically, based on a. linear variation of the 
parameter vector. The only non-zero inviscid source 
term is, 


B' s = P = 


7-1 


Z\Z\ 


z\ + z\ 


(16) 


The integration over the triangular element is divided 
into thirds along the median-dual boundaries, as in 
Figure 1, so that, 


fl — fli T O 2 T H 3 (17) 

The subintegrals are then distributed to the nearest 
node. Notice that the subdivided integration elements, 
fli _ 3 , are quadrilaterals, whereas the original element 
was a triangle. The distribution formula is thus, 

w aj Sj\Jj t f— w f Wdflj+COE (18) 

The integration of the source term over flj is expanded 
in detail in Ref. 11. 



Fig. 1 Subdivision of triangular element into three 
quadrilateral integration areas. Dashed lines are 
the median-dual mesh. 


The viscous axisymmetric source can be integrated 
using the Haselbacher 12 thin-layer approach (detailed 
in the following viscous flux section) as, 


/ 


bi <m 




2 

3i?eoo 



pV ■ hdT 

(19) 


Mass lumping to the node for the first term yields, 

f —dn =HiSi— (20) 

Jot y Vi 

while the second term is evaluated at edge midpoints. 

Inviscid Flux 

The finite volume discretization of the inviscid flux 
is performed as an average of the fluxes to the left and 
right of the control volume face, times the parameter 
w a (equal to 1 for two-dimensional or the y-value of 
the quadrature point on the face for axisymmetric) 
plus the artificial dissipation, which is defined as, 

* =i|l-n|(U R -U L ) (21) 


where by convention the right state is to the outside of 
the control volume while the left state is to the inside. 

The parameter vector, Z = [1, u, v, H] T , is 

linearly averaged, Z = |(Zl + Zr), to provide the 
quantities, 



and the Roe-density is, p = \fp\IpR- 

The projected flux Jacobian is decomposed as, 


|A-n| = X|A|X -1 (23) 

with. 


X = 


A 

= diag (V, 

V,V + a, 

V- 

a) 

1 

0 

1 


1 

u 

— n y 

u + an x 

u ■ 

- an x 

V 

n x 

v + an y 

v ■ 

- an y 

V 2 
, 2 

vn x — un y 

H + aV 

H 

— aV_ 


(24) 

(25) 
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The product X 1 (Ur — U|_) is expressed, 

( 2a? dp — 2 dP \ 
2 a?(n x dv — n y du) 
dP + pa dV 
\ dP-pddV ) 


X _1 dU = 


2 d 2 


(26) 


where the projected velocity is V = V -n and the aver- 
aged speed of sound for a perfect gas is, 


d 2 =(7-1) (H 


u 2 + v 2 


(27) 


Integration of the inviscid flux for fluctuation split- 
ting is performed as, 

f V-(rt 7 a F i ) dfl= [ WaV-^dfi + w [ F * d£l 
Jn ' ' Jn Jn 

(28) 

The y-component of the flux function can be written 
in terms of the parameter vector as, 


F* = 


(29) 


Z\ Z r i 

Z’lZ^ 

zi + 2=1 (ZiZ 4 - 

Z 3 Z 4 . 

A linear variation of the parameter vector over a tri- 
angular element can be represented as, 


Z{x,y) — 2 g^ijkZj 


(30) 


• [(« - Xi)(y k - yi ) + (y- Vi)(xi - x k )] 

where ey* is the cyclic-permutation summation oper- 
ator. The linear variation can also be written in the 
element-local (£, rj) coordinates, referring to Figure 2, 
as, 

Z (€,»?) = Z, + ^(Z 2 - ZOC + ^(Z 3 - Z 2 )», 

= Zi + — A^Z£ + —AjjZt] 

1 3 ti 

The domain is on 0 < ij < and 0 < £ < £ 3 . 
cartesian coordinates map similarly, 

x(€,v) =xi + ~A i x^+^-A v xr] 

13 1 1 

S/(£,»?) = 2/1 + 4 + T A nV V 

t-3 <4 

Some general integration rules can be developed for 
linear variations over the triangular elements: 


/ 

Jn 


f 

Jn 


dn =s T 

xdfl = Sjx 


£ xy d?l = .S T :cy - ^ ^;cy “ ! T x J y ^j 


(34) 

(35) 

(36) 



Fig. 2 Elemental triangular domain for fluctuation 
splitting. 


The cell-averaged value is, 


x = 


X\ + x 2 + X 3 _ 1 

“3 ~3^ X} 


= ,£** ( 37 ) 

3 = 1 


The last term of Eqn. 28 is distributed to the nodes 
in a manner similar to the source term, 

m ai SiUi t < zu ( F* y dfl, + COE (38) 

This integral can be evaluated exactly as, 


L z ' z > dS> ‘*m 


I 4 Z 1 Z 3 + 11(^1, Z 8 + Z\Z 3i ) 

(39) 


+ 9Zl; Z 3i + ^2 Zlj Z'ij 
i=i 

for the continuity equation. The integrals for tne other 
governing equations follow directly from Eqn. 39. 

The remaining term to evaluate in Eqn. 28 is the 
inviscid fluctuation, 


(31) 

0 = - [ WaV-F* dSl 




Jn 



The 

= -| 

) dfi 

(40) 


if/ 

\ 


(32) 

= J w o, [iihi-YzAifL - 

where dF* = Fz dZ and, 

l 3 n 3 -F zA n Z 

) dfl 

(33) 




F y z = 


2=1 Z 4 
7 4 

0 

0 

z?, 

0 


Zi 0 0 

2 ±iz 2 -J=zz 3 2=1 Zl 

-1/ xj y J- 


Z 3 

Z 4 

0 

Z 3 


7 

Z 2 

0 

Zr 

Z 2 


0 

Z 2 

0 

0 


(41) 


^4 -lfZ 2 2±I 

0 0 z 4 


^Z, 

Z 3 


(42) 
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The integration rule Eqn. 36 allows for the direct eval- 
uation of Eqn. 40 as, 


Similarity transformations lead to, 


(f) — CE7 0 ^ini- 


Fz-\ ( Fz-lE^-F Zj 
j=i 


A f Z 


+ 


3 = 1 


Noting that A = Fz and, 


A^Z 


Z/v = 


2VP 


1 

— u 
—v 


(43) 


0 0 
0 0 
2 0 


i?+(7 — 1) (« 2 +n 2 ) — 2(7 — l)w — 2(7 — l)v 27J 

(44) 

with the tilde-averaged quantities defined as, 

U = U(Z), A^U = U^A C Z, VJ = U 2 A,Z (45) 
leads to the fluctuation expressed as, 


(f) — ' 


+ 2 ro «^n3' 


A-jlA-'EfA, 

J =1 


A-ifA-iE^A, 

J =1 


= + (t> V 

Equation 46 includes the approximation, 


A i - Fz^o 


A f U 


A.,,U 


(46) 


(47) 


A transformation of variables can be made to the 
auxiliary variables so that, 


with, 


U w = 


<t> = U w<P 


100 4 

a 2 

u 1 0 \ -u 

a 2 

v 0 1 \v 

a 2 
V 2 


(48) 


u v 


-d— In , 

7-1 T J 


(49) 


0 — 2^ a^l^l* 




^ a 


4 y 3 jr[ 

4 l 3 ~i 


a A{W + (3 A,jWj 

= # + 4> n 

with A. = W(/AUff and, 

/Q\V = W [l A^J = W z A^Z 
A^W = Wj > A^U = W ? A„Z 


where, 


W 


1 / 


Wz = Vp 


o To 
Z T 

—U 

—V 

7rly2 
L 2 v 


2-4 

7 ft 

—u 

—v 


Z AA v 1 


A f W 


A^W 


(50) 

(51) 

(52) 


0 

■(7- l)w 


7/1 

1 


-(7-l)u 7-1J 


(53) 


11 

7/l 

0 

1 




7 


1 

7/1 

0 

0 

2=1 

7 ■ 


7 a 2 2 

Zi 

0 


70 2 A 

0 

Zi 


-^Zo -^Z, 


2Zi - 24 Z 4 

1 7 a 2 ^ 

-z 2 

-Zi 

^Z A 


using the relation, 

h =c p T= — - 

7-1 

The projection of .4. has the simple form, 

0 0 


-*=£ Zi 

7 a 2 1 

0 

0 

4^1 . 

(54) 

(55) 


n-A. = 


V 0 
0 V 
0 0 


0 

V 


0 a 2 n x a 2 n y V 


(56) 


where the projected velocity is V = n-V. The gener- 
alized advection speeds are, 


a =- ni - 


a 6 - 

P =-yn 3 - 


4 y 3 jr[ ®a 

4 y 3 


(57) 

(58) 
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where the approximation A.j = W^AjU^ has been 
incorporated. 

Linearity preservation for second-order spatial ac- 
curacy is obtained by limiting the fluctuations compo- 
nentwise, 

= $ + # i’iQj) = $ (i - 4' ) (59) 

4>f =4>]-4> n j 4’(Q j ) = 4>](i-4’(Q j )) (60) 

where the second equalities hold for symmetric lim- 
iters. The limiting ratio is, 


Qj 



(61) 


In vector form, Eqns. 59 and 60 can be written, 


with 


0* f =D«^ (62) 

**1 »_> T) 

<j> = D n <p (63) 

D ' =diag(l-^(2-)) (64) 

D" =diag(l -tp(Qj)) (65) 


Upwinding is achieved through the introduction of 
artificial dissipation 

4> = sign(a)^> = A^W (66) 

4>' V =sign(/3)0*” = -w a MgWMp 1 \/3\AfW (67) 

where M a = sign(a) and M.g = sign(/3). The ab- 
solute values of the generalized advection speeds are 
developed using the following decomposition, which is 
exact for the two-dimensional equations but approxi- 
mate for the axisymmetric form, 

a = ^fn-A= ^XAX* 1 (68) 

where A remains as in Eqn. 24 and, 


which for |V a | > a leads to, 

M a = sign (V a )J (73) 


0 0 
—n*n\ sign(V a ) % 
nf sign(V a ) 

an\ 0 

(74) 

where V a = h\-V. Similarly, defining |/3| as, 

m =-|<V|A|A” 1 (75) 

for | Vg 1 > a leads to, 

M 0 =sign (%)I (76) 

and for | | < a, 

0 o’ 

-n^n^sign(V /3 ) ^ 
nfsign(Vg) ^ 
anl 0 

(77) 


M« = 


sign(V/ 3 ) 0 

0 nfsign(Vg) 

0 -ngn^sign(V^) 


0 


am 


and for |V a | < d, 


= 


sign(V a ) 

0 

0 

0 


n\ sign(V a ) 
-n^sign(V a ) 


am 


where Vg = m-V. M a and have the property, 

M- 1 = M tt M/ = Mg (78) 

Eigen- value limiting for the suppression of expansion 
shocks can be introduced into Eqn. 70. If the limited 
eigen- value is expressed as, 

|A|, i TO HA| + <5 a (79) 



1 

0 

0 

0 

X = 

0 

n y 

n x 

h x 

0 

—n x 

n y 

h y 


_0 

0 

a 

—a 


The absolute value is then defined as, 

H = |^|A|A _1 (70) 

Expressions for the sign of the generalized 

wavespeeds are developed from, 

|a| = M a a (71) 

M a =A|A|A _1 A _1 (72) 


then the additional artificial dissipation for eigen- value 
limiting in the £ direction to be added to Eqn. 66 is, 
with 5 + = 1(<5a 3 + <5a 4 ) and d~ = |(<5a 3 - <5a 4 ), 

_ h 
Wa 2 ' 

6 Al 0 0 0 

0 nf <5 + + n\ 5\ 2 ^i«-i(<5 + — ^a 2 ) 

0 n^nl[((5 + — <5 a 2 ) n f <5a 2 + if <5 + ^-<5“ 

0 an\8~ an[5~ (5 + 

• aTiE (80) 

while the eigen-value limiting in the t] direction takes 
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the form, 


analytically for a linear variation of the parameter vec- 
tor, 


, - 4 

+ H7 a — 


S\ 1 

0 

0 

0 

0 

rig 2 (5 + + rig <5 a 2 

n x 3 n y 3 (S + -5 a 2 ) 


0 

n§n^(<5+ -Sa 2 ) 

nf 5 a 2 + n \ <5 + 

a 

0 

an 3 S~ 

an y 3 S~ 

6+ 


•A V W (81) 

The fluctuation is distributed to the nodes using, 


3 

= - ^ E ^ = ~ £ ^ ( 88 ) 

3 = 1 

<t>i =^rJ Q -fiidn = ^ (89) 

Struijs et a/ 14 have shown that derivatives of primi- 
tive variables can be consistently represented in terms 
of the parameter vector gradients as, 


SiUu 

•*- 

-4*) + COE 


s 2 u 2t 

•*- 


+ </>'") + COE 

s 3 u 3t 

•*- 

-<f) + COE 

(82) 


where COE stands for contributions from other ele- 
ments. The distribution can be expressed in a more 
compact form, 


SjU* 


i(3-*)(0* 4 +(-!)>'*) 


+ (-4 + 5f-* 2 )(0*’’ + (-l) i 0'”) 
+ COE * = 1,2,3 


(83) 


where Eqn. 48 is used to define, 

f 5 =%0 t£ , <t>' ( =\Jw4>' i (84) 

<t>* v =Vw4>*\ </>'"= Ui v4>' V (85) 


Viscous Flux 

Viscous terms are discretized using a non-upwind 
fluctuation splitting formulation, which is equivalent 
to a finite element discretization using mass-lumping 
to the nodes. Both the finite volume and the upwind 
fluctuation splitting inviscid discretizations have been 
shown by the authors 13 to be compatible with this 
viscous treatment for scalar equations. 

Integrating the viscous flux over a triangular ele- 
ment leads to the viscous fluctuation, 


VV =V^VZ 


(90) 


where. 


Vz = 


and, 


2Zi 

0 

0 

0 

Aa 

i 

Zl 

0 

0 

"ft 

0 

l 

Zi 

0 

. 7 4 

7 z 

.m 

N 

74 

H 

1 

3^1 Zi 

7 1 


v =V(Z) = 


^l 2 

la 

Zi 

la 

Zi 


V 2 ^ [ Z ^ Z ^ ~h { Z 2 + Z l)\) 


(91) 


(92) 


Further, the consistent temperature gradient is, 


VT = 


VP PVp _ 
pr p 2 r ~ up 2 y 


(pVP-PVp) (93) 


The viscous fluctuation is then distributed to the 
nodes, 

SiV it <- tf + COE (94) 

An alternate approach to integrating the viscous 
flux is obtained by using the divergence theorem, 



= WaF v -hdY 


(95) 


<f) v = f V-{w a F v )dfl (86) 

Jn 

The nodal distributions are developed in a finite ele- 
ment sense by integrating by parts, 

<f>i = <f Viw a F v hdT — [ w a F v ■ dfl (87) 
Jr Jn 

For interior nodes the boundary integral in Eqn. 87 
will sum to zero and the volume integral is integrated 


where fl, is the generalized control volume containing 
node i, with two-dimensional area equal to Si, and T,; 
is the boundary of f V 

Haselbacher et al 12 have recently presented an ap- 
proximate treatment for integrating Eqn. 95 on two-di- 
mensional unstructured grids, which they relate to the 
thin-layer approximation of the Navier-Stokes equa- 
tions presented by Gnoffo 15 for structured grids. The 
method preserves positivity for the Laplacian and is 
transparent to grid topology. 
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The development begins with the expression, 


¥°-h = 


R e 


Vw-h + |V-Un x 


Vv-t 


Vv-n+ |V-Un y + Vu-t 


KVT-h + M |§(V-U)(U-n) 
y+^Vy 2 -h — uWv-t + vVti-i 


(96) 


/ 


where t is a tangent vector with components (— n y , n x ). 
Haselbacher’s approximation neglects all tangential 
terms and approximates V ■ V ~ n x Vu-h + n y Vw-n. 
Using the notation it n = Vie-n, etc., and including the 
axisymmetric terms gives the approximation, 


leading to, 


V-U ~ Vn'fl + G7- 

y 


(97) 


F”-n 



( ° \ 

K T + W ( V n -h + w 

V ^T„ + U-U„ + |U-h (U n -h + w v ~) j 

(98) 



Interior 


Exterior | 

■ =Boundary state, 


= Reconstructed state 


Fig. 3 Weak implementation of finite volume 
boundary condition for node 0, imposed by spec- 
ifying external state. Quadrature points denoted 
by X’s. 



Fig. 4 Weak implementation of fluctuation split- 
ting boundary condition, imposed by specifying 
external state at ghost nodes fo, f\. 


A further simplification aligns n with the nearest 
mesh edge for faces of T on the interior of the domain, 
so that terms such as u n reduce simply to the difference 
in nodal values divided by edge length. Also, w a is 
chosen to be the midpoint of the edge. 

Including one more approximation, namely replac- 
ing the length T of the median-dual face with the 
length of the associated containment-dual face, has 
the effect of canceling some of the errors for very 
high-aspect-ratio cells introduced by assuming n is 
edge-aligned. For low-aspect-ratio cells, the contain- 
ment dual is the same as the median dual and the true 
n is closely aligned with the mesh edge. This imple- 
mentation is similar to the suggestions of both Barth 6 
and Haselbacher, 16 yet retains the global median-dual 
implementation required by a distribution scheme. 

Boundary Conditions 

Boundary nodes may be updated either strongly, 
where the nodal solution values are simply assigned 
a priori , or weakly, where the solution values at the 
boundary nodes are relaxed in time using the same 
formulations as for interior nodes. 

For finite volume, the weak boundary implemen- 
tation specifies the solution state to the outside of 
boundary faces, then allows the approximate Riemann 
solver to construct the appropriate fluxes through the 
boundary faces. See Figure 3 for an illustration of the 
weak finite volume boundary condition. The solution 
state to the inside of the boundary face can be set from 
either a first- or second-order reconstruction from the 


node. For some cases, second-order reconstruction to 
boundary faces has led to localized oscillations in the 
solution convergence at boundary nodes. 

Weak formulation of the fluctuation splitting bound- 
ary condition is developed using fictitious “ghost” 
nodes, one for each boundary node, as shown in Fig- 
ure 4. Considering the scalar case, the positioning of 
a ghost node such that the edge connecting the ghost 
and boundary nodes is parallel to the advection veloc- 
ity results in a boundary fluctuation of, 

fibco = 2^ 01 ^'^ 01 (^/o — ^o) (99) 

for node 0 of Figure 4. Observe that this formulation is 
independent of the physical location of the ghost node, 
so the ghost node can be chosen to be infinitesimally 
close to the boundary node it supports. The solution 
state at the ghost node remains to be specified, and 
can be varied node-to-node. The associated artificial 
dissipation is, 

fibco = sign(Xn 0 i)^6c 0 (100) 

and the resulting distribution is, 

SoU()t -ji^bco + 4>bc 0 ) + COE (101) 

Since no account of the ghost cell area is made in form- 
ing the dual area on the LHS of Eqn. 101, a scale factor 
on [|, 1] can be applied to the distribution. 
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The extension to systems follows by analogy. The 
boundary fluctuation is defined, 


Both the full viscous flux and the approximate thin- 
layer flux of Haselbacher reduce to, 


he o = ^oiAc 0 -noi(W /o - Wo) (102) 


F v -n = 


R e 


0, P 1 n . K,T n 


(107) 


with the artificial dissipation. 


4>bc n = sign(^ C0 -hoi )<t> bco 


(103) 


at a wall, since V, V-V, and V n -n go to zero. Defining 
the heat transfer into the wall according to Fourier’s 
law, 


Additional dissipation for the suppression of expan- 
sion shocks is added to Eqn. 103 following the form of 
Eqns. 80 and 81 as, 


Qw 


K 

r7 c 


T 

n 


and the wall shear stress, 


(108) 
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•(w /o 

-Wo) 


(104) 


The distribution to the boundary node is then formed 
as, 


5 0 U 0t «-^U- w(i beo +& eo ) + COE (105) 

This system treatment is only approximate, as the 
cross-fluctuation does not vanish when V || / o 0, as in 
the scalar case, but reduces to the term, 


0 0 0 0 

0 0 0 n% 

0 0 0 n y fo 

0 o 2 «-/ 0 a 2 n^ o 0 


(106) 


Strictly, there should be some cross-coupling with the 
neighboring boundary nodes. However including the 
term from Eqn. 106 requires explicitly locating the 
ghost nodes, which can be impossible for certain ge- 
ometries. 

The distribution to node 1 is formed analogously, 
substituting 0 for 1 in Eqns. 99-106. 

The freestream boundary condition is enacted by 
specifying a complete, constant thermodynamic state 
and velocity vector. By using the weak boundary en- 
forcement, this one boundary condition covers the four 
permutations of subsonic or supersonic, inflow or out- 
flow. 

The inviscid wall boundary is implemented by mir- 
roring the primitive variables, either across the face for 
finite volume or at the node for fluctuation splitting. 

The axisymmetric axis is enforced by imposing zero 
flux on the axis and using the control- volume centroid 
for w a in Eqn. 8. 

Viscous walls define a stagnant velocity and a spec- 
ified wall temperature. The zero velocity at the wall 
causes the viscous axisymmetric source to be zero. 


f — ^ V 

• w — D v n 

-^eoo 

(109) 

allows the Eqn. 107 to be written as, 


F * ft — [0, T w , Qw] 

(110) 


where the minus sign results from the choice of an 
outward unit normal, n, to the control volume, which 
points into the wall at a boundary. 

The solid wall is enforced weakly, by specifying the 
wall shear that will drive the flow momentum to zero 
and the heat flux that will drive the solution temper- 
ature to the desired wall temperature. An advantage 
of this weak approach is that wall heat transfer and 
skin friction are solved for directly, rather than as a 
post-processed least-squares reconstruction. Using an 
explicit update, the wall heat flux can be isolated as, 


Qw — 




ZUaS 

At 


(U 4 - L\ e(T w )) + RHS\ +V 


(111) 


Similarly, the wall shear is, 


Tin — 




ZU a S 

~KT 


(U 2 ,U 3 ) + RHSi+ v 


( 112 ) 


Temporal Evolution 

Solutions at the nodes are updated using an explicit 
forward-Euler LHS. A Jacobi relaxation strategy is 
followed with either local time stepping or first-order 
global time steps. 

The CFL (Courant et al 17 ) criteria for explicit 
schemes is adapted for use with the node-based un- 
structured scheme. The inviscid timestep is defined 
by the most restrictive time for signal propagation, at 
the eigen-value speeds, between adjacent nodes, 


Ato = min 




roi-roi 


Vl^o-roil + a o 


|Uo-roi| + ao ||fbj|| J 

(113) 


where the current node is node 0 and i takes on nodal 
values for each distance-one neighbor of the current 
node. 
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Y 


Fig. 5 Sample randomly distorted mesh used for 
solver verification cases. 



The viscous timestep restriction is taken to be an 
approximation based upon the positivity analysis for 
the scalar case, 13 assuming order-1 Prandtl number, 


Ato 


4:S 0 p 0 R e<xi R 
fi 0 (R + 7 - 1) St 


(114) 


The stability and convergence of axisymmetric solu- 
tions is found to be enhanced by scaling the timestep 
for points near the cylindrical axis by the maximum 
of either the node height or the square root of the 
median-dual area. 

The more restrictive of the inviscid or viscous 
timestep is used to scale the nodal update. 


Results 

Verification of the complete solver is performed in 
stages using a methodology derivative of Singhal. 1 A 
variety of canonical cases are constructed, including 
grid distortions, that are designed to exercise com- 
binations of the various functions that comprise the 
complete solver. Two viscous cases serve to validate 
the solver against benchmark data. 

Inviscid Verification 

Distorted mesh 

The first verification case simply passes a uniform 
flow through a distorted grid, with success being the 
preservation of uniformity to at least six significant 
digits. The domain is initialized to stagnant condi- 
tions with freestream flow impulsively applied at the 
boundaries. A variety of flow angles were tested on 
— 180° < AOA < 180° for subsonic, transonic, and 
supersonic Mach numbers. Regular, high aspect ratio 
(100), skewed (2° < 9 < 175°), and randomly dis- 
torted (Figure 5) meshes with 121 nodes were used. 
Initial runs were instrumental in refining the treat- 
ment of eigen-value limiting for fluctuation splitting. 
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State A 
Mach = 2.3 
AOA = -10 


State B 
Mach =1.8 
AOA = 10° 


State C 
Mach = 1.91 
AOA = 0° 



State D 
Mach = 1.45 
AOA = 0° 


Fig. 6 Description of converging-Mach-stream 
problem. Flow from left to right, with oblique 
shocks, solid, and slip-line, dashed, emanating from 
trailing edge of splitter plate. 


All final runs were successful for both finite volume 
and fluctuation splitting. 

Converging Mach streams 

Thermodynamic routines are verified by considering 
converging Mach streams, inclined at ±10°. The up- 
per stream is at Mach 2.3 while the lower stream has 
Mach 1.8. The two streams have matched densities 
but a temperature ratio of 1.0812, resulting in a hori- 
zontal slip line behind the oblique shocks. A complete 
description of the analytic solution appears in Figure 6 
and Table 1. 

A sequence of four meshes, with a refinement ratio 
of 1.5, is considered. The meshes are triangulated from 
16 x 16, 24 x 24, 36 x 36, and 54 x 54 grids. The tri- 
angulated 16 x 16 grid is shown in Figure 7. The finer 
meshes cover the same domain and are constructed 
similarly to the shown mesh. 

A Mach-number contour plot for fluctuation split- 
ting on the finest mesh is shown in Figure 8, showing 
crisp discontinuity resolution and the correct post- 
shock Mach numbers. The shock angles for all eight 
cases, i.e. finite volume and fluctuation splitting on 
each mesh, are measured to be correct within ±1°. 
The ^ 2 -norms of the primitive- variables error at states 
C and D are plotted versus the characteristic mesh 
size in Figure 9. The slopes of the regression lines 
are indicative of the order of accuracy with respect to 
grid convergence of the two schemes for this test case. 

Table 1 Analytic thermodynamic states for con- 
verging Mach streams. 


State 

P , kg/m 3 

T, IC 

P, kPa 

V, m y 

A 

1.2 

300 

103.34 

798.6 

B 

1.2 

324.7 

111.73 

649.9 

C 

1.813 

356.7 

185.62 

723.7 

D 

1.718 

376.4 

185.62 

563.7 
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10° 



Fig. 7 16 x 16 mesh for converging Mach streams. 



Fig. 8 Mach contours from fluctuation splitting 
on finest mesh. 


Fig. 10 Convergence histories for converging 
Mach-stream case. Fluctuation splitting solid, fi- 
nite volume dashed. Coarsest mesh on left, finest 
on right. 



Fig. 11 Grid for diamond airfoil verification test. 


logio \\Error\\ _ 


2:1 slope 


3:1 slope 


log 10 V nodes 


fluctuation splitting and 2.1 for finite volume. 

Temporal convergence rates are plotted in Figure 10, 
with timings performed on an IRIS R10000 platform. 
All cases were run using the Minmod limiter and a 
Jacobi update strategy with local time steps. Fluc- 
tuation splitting was run with a unity CFL number, 
while best convergence for finite volume was found for 
CFL=0.7. Fluctuation splitting runs at 145 /xs per 
node per iteration, while finite volume runs at 165 /xs 
per node per iteration. 


Fig. 9 Grid convergence rates for converging 
Mach stream case. Circles = fluctuation splitting, 
squares = finite volume. 

Finite volume exhibits second-order convergence, as 
expected. Unexpectedly, fluctuation splitting shows 
super-convergence for this particular case. True mul- 
ti-dimensional upwinding is likely the source of the ex- 
ceptional fluctuation splitting accuracy for this purely- 
supersonic flow. Supplementing the graphical deter- 
mination of the grid-convergence rates, the equations 
presented by Roache, 2 based on a Richardson extrap- 
olation, yield average grid-convergence rates of 3.0 for 


Diamond airfoil 

A verification of the inviscid wall boundary condi- 
tion is performed on a diamond airfoil at zero angle 
of attack and Mach 1.5. The flow deflection is five de- 
grees. The grid is shown in Figure 11. A Mach- number 
contour plot using fluctuation splitting is shown in Fig- 
ure 12. The corresponding finite volume solution, not 
shown, is visually indistinguishable from the fluctua- 
tion splitting solution. The analytic drag coefficient, 
based on chord length, is 0.02760. The fluctuation 
splitting drag coefficient is 0.02638, for a 4.4 percent 
error. The finite volume result has an error of 6.6 per- 
cent from a drag coefficient of 0.02579. 
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Fig. 12 Mach contours on diamond airfoil, M = 1.5, 
fluctuation splitting solution. 



Fig. 13 Two-dimensional 10 percent circular- 
bump mesh with isobars from fluctuation splitting 
solution at Mach 0.1. 


Circular bump 

A subsonic two-dimensional verification is per- 
formed on a 10 percent circular bump at Mach 0.1. 
The 1389-node mesh with isobars from the fluctuation 
splitting solution is shown in Figure 13. A true incom- 
pressible inviscid flow would have symmetric isobars 
fore and aft, and zero drag. The fluctuation splitting 
drag coefficient, based on cord length is 0.0058. Finite 
volume predicts a drag coefficient more than twice as 
large, 0.0128. A lower fluctuation splitting drag coeffi- 
cient is indicative of lower levels of artificial dissipation 
in the solution for this case. 

Sphere 

In a similar vein, Mach 0.1 axisymmetric flow over a 
sphere is tested on a 1369-node mesh. The drag coef- 
ficient, based on frontal area, is 0.43 for finite volume 
but 0.56 for fluctuation splitting. Contrary to expec- 
tation, the increased artificial dissipation in the finite 
volume solution creates enough of a total pressure loss 
to nearly eliminate separation on the leeside, whereas 
the leeside increase in pressure toward the centerline 
in the fluctuation splitting solution does produce a siz- 
able separation region, and in this case a larger drag 
coefficient. As with subsonic bump case, true incom- 
pressible, inviscid flow should theoretically produce 
zero drag. 
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Cone 

The final inviscid verification is for an 11-degree 
semi- vertex-angle cone at Mach 1.5. The well- 
established Taylor-Maccoll 18 method for the conical 
supersonic Euler equations predicts a drag coefficient, 
based on base area with no base pressure, of 0.7795. 
The fluctuation splitting solution, which converged 
seven orders of magnitude in 38 seconds, predicts a 
drag coefficient of 0.7785, for only a 0.13 percent er- 
ror. The finite volume solution, which took 13 percent 
longer at 43 seconds to reach seven orders of magnitude 
residual convergence, predicts a 0.7754 coefficient, for 
an error of 0.53 percent. 

Viscous Validation 

Two canonical viscous validation cases are consid- 
ered: a subsonic flat plate and a hypersonic cylinder. 
Steady laminar solutions are obtained using the Hasel- 
bacher thin-layer viscous treatment with containment- 
dual modification. 

Flat plate 

The classic Blasius 19 flat-plate boundary layer prob- 
lem is solved on a rectangular domain. Mach 0.3 
flow enters 2 units upstream of the plate leading edge, 
which is located at the origin. The plate is 4 units 
long, ending at an extrapolation outflow boundary. 
The upper boundary is 1.2 units above the plate. The 
Reynolds number is 10 4 . 

The meshes are obtained from a structured grid con- 
taining 37 equally-spaced points parallel to the plate, 
12 points upstream of the plate and 25 points on the 
plate, and 41 points normal to the plate. The ver- 
tical grid spacing at the wall is 0.001 units with an 
exponential stretching as described in Ref. 20, plac- 
ing approximately 20 nodes within the boundary layer. 
The unstructured mesh is formed from the structured 
grid using diagonal cuts in an alternating pattern. 
Two coarser meshes are similarly constructed by suc- 
cessively deleting every-other node in the wall-normal 
direction, leaving 10 and 5 nodes, respectively, in the 
boundary layer for the medium and coarse grids. 

Boundary layer profiles of u are extracted at x = 
1,2,3 from both the fluctuation splitting and finite 
volume solutions and plotted versus the Blasius solu- 
tion in Figure 14. The boundary layer scaling variable 
is defined as, 



Both solution sets match the Blasius profile, indicating 
well-developed flow with adequate grid resolution on 
the finest mesh. 

Figure 15 shows the effect of using the containment- 
dual approximation in the Haselbacher thin-layer vis- 
cous treatment. Boundary layer profiles of u are again 
extracted at x = 1,2,3, with both solutions being run 
with fluctuation splitting. Figure 15(a) is the same 
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a) Fluctuation splitting. b) Finite volume. 

Fig. 14 Boundary layer profiles of tangential velocity extracted from three stations on flat plate. 



a) Containment dual approximation. b) Strict median dual implementation. 

Fig. 15 Boundary layer profiles computed using two different viscous dual mesh definitions. 
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Fig. 16 Boundary layer profiles of vertical velocity 
extracted from midpoint of flat plate. 


as Figure 14(a), while Figure 15(b) uses the strict 
median-dual definition for the viscous terms. For the 
highly-stretched grid elements used in this case, it is 
clear that the containment-dual approximation pro- 
vides improved boundary-layer resolution, while omit- 
ting the approximation leads to a profile that is “too 
full.” 

The u-velocity profiles from the fluctuation splitting 
and finite volume solutions are compared in Figure 16, 
both extracted from the plate at x = 2. The fluctua- 
tion splitting solution comes much closer to matching 
the Blasius profile than the finite volume result. Ex- 
cessive artificial dissipation is produced by the finite 
volume scheme in the .(/-momentum equation, which 
suppresses the u-velocity below the analytic value. The 
artificial dissipation contributions to the y-momentum 
equation are plotted for both fluctuation splitting and 
finite volume in Figure 17. The vertical scale has been 
enlarged by a factor of 30 to zoom in on the boundary 
layer in Figure 17. Clearly, finite volume is producing 
significantly more artificial dissipation than fluctua- 
tion splitting over the length of the boundary layer. 

For this essentially incompressible case, the suppres- 
sion of the vertical velocity due to excessive artificial 
dissipation is manifested by an increase in skin friction 
coefficient, as shown in Figure 18, where the friction 
coefficient increases with running length for finite vol- 
ume, but not for fluctuation splitting. Recall that 
finite volume is continuously producing artificial dissi- 
pation over the length of the plate while the fluctuation 
splitting dissipation is restricted to the leading-edge 
region only. Figure 18 presents data from all three 



b) Finite volume. 

Fig. 17 Artificial dissipation production in the 
(/-momentum equation. Eleven contours spaced 
equally on 0-0.0005. 


grid refinement levels. The finite volume results de- 
grade dramatically with coarsening of the mesh, but 
the fluctuation splitting results remain relatively in- 
variant with mesh resolution, all the way down to only 
five nodes in the boundary layer. 

The medium-mesh finite volume solution was re- 
peated using the full Navier-Stokes treatment, rather 
than the thin-layer equations. No change in the skin- 
friction results are seen over the first half of the plate, 
Figure 19, though there is an 8-percent improvement 
toward the end of the plate. Solving for the full Navier- 
Stokes terms requires 11 percent more CPU time per 
iteration. 

Cylinder 

The opposite end of the Mach-number spectrum 
is used to validate heat-transfer calculations, in this 
case for a cylinder of 1 m radius in Mach 17.6 flow. 
The perfect-gas assumption is a poor physical model 
for these extreme conditions, V/o = 5 krn/s, = 
0.001 kg/m 3 , Too = 200 K, T wa u = 500 K, but the 
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b) Finite volume. 

Fig. 18 Skin friction coefficients for Blasius flow. 
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tered to the wall, and circumferential nodes spaced 
every 3 degrees. Only the forward-half of the cylinder 
is solved, as shown in the mesh and flowfield solution 
of Figure 20. 

The surface pressure coefficient is plotted versus ro- 
tation angle from the stagnation point for both the 
fluctuation splitting and finite volume solutions, along 
with the LAURA and FUN2D results in Figure 21. 
The LAURA, FUN2D, and fluctuation splitting curves 
all over-plot, and the finite volume solution nearly 
over-plots, being 1 percent low at the stagnation point 
and slightly high by a similar amount 90 degrees away. 
The calculations were repeated on a grid coarsened by 
a factor of four (skip of two in both structured-grid 
directions), with surface pressure results plotted in 
Figure 22 along with the fine-mesh LAURA solution. 
The coarsened fluctuation splitting surface pressures 
retain good agreement, and the finite volume solution 
matches over most of the cylinder, with minor excep- 
tions again at the stagnation point, 1 percent high on 
this grid, and at the 90 degree point. 

Surface heat-transfer rates for LAURA, FUN2D, 
and fluctuation splitting are shown in Figure 23. Both 
of the unstructured-mesh solutions show elevated heat- 
ing at the stagnation region, with fluctuation splitting 
being 30 percent higher than LAURA while FUN2D 
is 50 percent higher. Heating results for this case were 
also obtained using a validated structured-mesh cou- 
pled inviscid/boundary-layer code, with the solution 
agreeing with the LAURA data. 

The fine-mesh solutions were repeated using the full 
Navier-Stokes treatment, and no changes in heating 
levels were observed. 
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Fig. 19 Effect of viscous modeling on skin friction. 


case provides a severe test of the algorithms under a 
re-entry scenario. Results are compared against the 
LAURA 15 ’ 21 ’ 22 benchmarks.* The LAURA code is 
well-established as a structured-grid hypersonic solver. 
Also included in the LAURA benchmark data is a solu- 
tion using the unstructured-mesh finite volume solver 
FUN2D. 23 The unstructured grid for this case was 
obtained by simple triangulation of the LAURA grid, 
which has 65 nodes perpendicular to the surface, clus- 

*http: //hef ss . larc.nasa.gov 


Concluding Remarks 

A multi-dimensional upwind fluctuation splitting 
scheme has been formulated for two-dimension- 
al/axisymmetric viscous flows. A weak form of the 
boundary conditions was proposed and the proper in- 
corporation of eigen-value limiting derived. 

While the scheme is formally second-order accu- 
rate, super-convergent third-order behavior was seen 
for a canonical verification test of converging super- 
sonic streams. 

The fluctuation splitting scheme produced more ac- 
curate solutions for the inviscid diamond airfoil and 
circular bump than a finite volume scheme. Para- 
doxically, the excessive artificial dissipation produced 
by the finite volume scheme actually led to a lower 
drag than fluctuation splitting for the subsonic invis- 
cid sphere case. 

For the viscous flat plate, fluctuation splitting was 
seen to produce more accurate solutions than finite 
volume, due to the fluctuation splitting low levels of ar- 
tificial dissipation. Also, fluctuation splitting showed 
excellent skin-friction predictions on extremely coarse 
meshes, while the finite volume results deteriorated 
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0 30 60 90 

Rotation degrees from stagnation point 

Fig. 21 Cylinder surface pressures, solid = fluctu- 
ation splitting, LAURA, and FUN2D, while dashed 
= finite volume. 



0 30 60 90 

Rotation degrees from stagnation point 

Fig. 22 Cylinder surface pressures on coarsened 
mesh, solid = LAURA, dashed = fluctuation split- 
ting, and dotted = finite volume. 



b) Pressure contours. 

Fig. 20 Hypersonic cylinder domain with fluctua- 
tion splitting solution. 


Fig. 23 Cylinder surface heat-transfer rates, 
solid = LAURA, dashed = FUN2D, and dotted 
= fluctuation splitting. 
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with mesh coarsening. 

Surface pressures were well predicted for the hyper- 
sonic viscous cylinder, but surface heating was dis- 
appointingly high for all the unstructured schemes 
considered here relative to the benchmark solution. 
Further testing is warranted to probe the validity 
of unstructured meshes for hypersonic viscous so- 
lutions. Mixed-element strategies may be a more- 
appropriate course for high-Mach-number applica- 
tions, with the fluctuation splitting inviscid distribu- 
tion assigned from an implicit triangulation in the 
boundary layer. 

Therefore, the upwind fluctuation splitting invis- 
cid discretization is an attractive solver from subsonic 
through hypersonic regimes vis a vis finite volume for a 
fresh code build, but the benefits are perhaps not great 
enough to justify rebuilding a legacy code. However, 
when coupled with a viscous discretization the reduced 
levels of artificial dissipation in fluctuation splitting 
allow for coarser viscous meshes, which can lead to 
significantly reduced computational requirements. 
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